Analysing spatial patterns of the landscape
Intro2R-spatial workshop
author: Song, Xiao Ping
updated: 2020-03-23
1 Set up
Install necessary packages:
install.packages("tidyverse", dependencies = TRUE)
install.packages("raster", dependencies = TRUE)
install.packages("sf", dependencies = TRUE)
install.packages("landscapetools", dependencies = TRUE)
install.packages("landscapemetrics", dependencies = TRUE)2 Prepare data
2.1 Quick visualisation
Import only the RGB color bands as individual RasterLayer objects:
library(raster)
#blue
b2 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B2.tif')
#green
b3 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B3.tif')
#red
b4 <- raster('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B4.tif')Combine the RasterLayer objects and visualise the satellite image:
landsatRGB <- stack(b4, b3, b2) #order is impt
plotRGB(landsatRGB,
stretch = "lin") #scale the values (try using "hist" also)
2.2 Import data
Import all 5 bands from the satellite data as a RasterStack object named landsat:
filenames <- paste0('data/Landsat 8 OLI_TIRS C1 Level-1/LC08_L1TP_125059_20180524_20180605_01_T1/LC08_L1TP_125059_20180524_20180605_01_T1_B', 1:5, ".tif")
landsat <- stack(filenames)
#rename bands
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR')Check coordinate reference system of landsat:
crs(landsat)## CRS arguments:
## +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,02.3 Crop data
Import polygon of city boundaries as sgshp and check if their coordinate reference systems match:
library(sf)
sgshp <- st_read("data/master-plan-2014-region-boundary-web-shp/MP14_REGION_WEB_PL.shp")Check coordinate reference system of sgshp:
crs(sgshp)## [1] "+proj=tmerc +lat_0=1.366666666666667 +lon_0=103.8333333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs"Transform sgshp to the match the coordinate reference system of the landsat:
sgshp <- st_transform(sgshp, crs = crs(landsat))Crop landsat according to city boundaries:
landsat <- crop(landsat, sgshp) #crop to rectangle
landsat <- mask(landsat, sgshp) #mask values according to shape of sgshpPlot the cropped image using only the RGB bands:
landsatRGB <- subset(landsat, c(4,3,2)) #Red, Green, Blue
plotRGB(landsatRGB,
stretch = "lin")
3 Classify land cover
3.1 Calculate NDVI
Create a function that calcuates the Normalized Difference Vegetation Index (NDVI) for each pixel:
ndvi <- function(x, y) {
(x - y) / (x + y)
}Apply function to the NIR and Red bands of landsat
landsatNDVI <- overlay(landsat[[5]], landsat[[4]],
fun = ndvi)Limit the range of values to be from -1 to 1:
landsatNDVI <- reclassify(landsatNDVI, c(-Inf, -1, -1)) # <-1 becomes -1
landsatNDVI <- reclassify(landsatNDVI, c(1, Inf, 1)) # >1 becomes 13.2 Visualise NDVI
Map out the NDVI values:
plot(landsatNDVI,
col = rev(terrain.colors(10)),
main = "Landsat 8 NDVI",
axes = FALSE, box = FALSE)
Plot histogram of NDVI values:
hist(landsatNDVI,
main = "Distribution of NDVI values", xlab = "NDVI",
xlim = c(-1, 1), breaks = 100, yaxt = 'n')
abline(v=0.2, col="red", lwd=2)
abline(v=0, col="red", lwd=2)
3.3 Classify using NDVI
Create a matrix to be used as an argument in the reclassify() function:
reclass_m <- matrix(c(-Inf, 0, 1, #water
0, 0.2, 2, #urban
0.2, Inf, 3), #veg
ncol = 3, byrow = TRUE)
reclass_m## [,1] [,2] [,3]
## [1,] -Inf 0.0 1
## [2,] 0.0 0.2 2
## [3,] 0.2 Inf 3Classify land cover using the defined threshold values:
landsatCover <- reclassify(landsatNDVI, reclass_m)Plot the land cover classes:
3.4 Save raster
Save the reclassified raster landsatcover in the GeoTiff format:
writeRaster(landsatCover,
filename = "clean_data/landsat_landcover.tif",
overwrite = TRUE)4 Landscape metrics
4.1 Quick visualisation
library(landscapemetrics)
library(landscapetools)
#landsatCover <- raster('clean/landsat_landcover.tif') #reload saved raster
show_landscape(landsatCover, discrete = TRUE)
Check if the raster data is in the right format:
check_landscape(landsatCover)## layer crs units class n_classes OK
## 1 1 projected m integer 3 ✓4.2 Patch-level
E.g. Area of each patch in the landscape:
lsm_p_area(landsatCover)## # A tibble: 12,580 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 patch 1 1 area 14.6
## 2 1 patch 1 2 area 0.09
## 3 1 patch 1 3 area 27.3
## 4 1 patch 1 4 area 0.09
## 5 1 patch 1 5 area 0.18
## 6 1 patch 1 6 area 0.36
## 7 1 patch 1 7 area 0.09
## 8 1 patch 1 8 area 0.09
## 9 1 patch 1 9 area 0.09
## 10 1 patch 1 10 area 0.09
## # … with 12,570 more rows4.3 Class-level
E.g. For each class, the total area of all patches:
lsm_c_ca(landsatCover)## # A tibble: 3 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA ca 6742.
## 2 1 class 2 NA ca 28496.
## 3 1 class 3 NA ca 42933.E.g. For each class, the average area of patches:
lsm_c_area_mn(landsatCover)## # A tibble: 3 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA area_mn 3.22
## 2 1 class 2 NA area_mn 4.52
## 3 1 class 3 NA area_mn 10.34.4 Landscape-level
E.g. Total area of the landscape (all three land cover classes):
lsm_l_ta(landsatCover)## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA ta 78170.5 Credits
Spatial data used in this document:
- Landsat-8 satellite images from the U.S. Geological Survey
- Singapore Regional Master Plan 2014 from the Urban Redevelopment Authority
Copyright (c) 2020 Song, Xiao Ping